################################################################################
## Title   : Bases eleitorais da Segurança Pública em 2018 
## Subtitle: violência e capacidades locais
##
## Created: 2023-07-13
## Updated: 2024-04-24
################################################################################

####################################################
## 01. Set Directory Structure                    ## 
####################################################

getwd()   
path <- "D:/2024_bpsr_candidatos_da_segurança" # Change to your working directory 
setwd(path)
getwd() 

####################################################
## 02. Set Environment                            ## 
####################################################

library("pacman")
pacman::p_load(readxl,data.table,                            
               tidyverse, dplyr, skimr,stargazer, labelled,   
               writexl,corrplot, margins,sandwich,knitr,fixest,                                    
               car, ggcorrplot, coefplot, texreg, prediction, 
               lmtest, sandwich, miceadds,margins, interplot, 
               sf, sp, spdep,geobr, grepel, gridExtra)

####################################################
## 03. Loading external data                      ## 
####################################################

dados <- read_excel("./dados/data_bpsr_analysis.xls")
glimpse(dados)  # Visualize the dataframe structure
colnames(dados) # Visualize the dataframe structure (columns names)

####################################################
## 04. OLS models                                 ## 
##                                                ## 
## Article: Table 1                               ##
####################################################

models <- list()

models$m1 <- lm(VOTOS_SEGURANCA_PROPORCAO ~
                  HOMICIDIO_2018_LOG             + factor(CAPACIDADE)+
                  CANDIDATOS_SEGURANCA_PROPORCAO + NEP               +
                  MAGNITUDE                      +
                  PREFEITOS_ESQUERDA             + PREFEITOS_CENTRO  +                  
                  ESQUERDA_2014                  + CENTRO_2014       + 
                  PIB_PERCAPITA_LOG              + POPULATION_LOG    +
                  URBANIZACAO                    + EVANGELICOS       +
                  NAOBRANCOS_2010                + POBRES            + 
                  GINI                           + 
                  factor(REGIAO), 
                data=dados)

models$m2 <- lm(VOTOS_SEGURANCA_PROPORCAO ~
                  HOMICIDIO_2018_LOG             + factor(CAPACIDADE)+
                  CANDIDATOS_SEGURANCA_PROPORCAO + NEP               +
                  MAGNITUDE                      +
                  PREFEITOS_ESQUERDA             + PREFEITOS_CENTRO  +                  
                  ESQUERDA_2014                  + CENTRO_2014       + 
                  PIB_PERCAPITA_LOG              + POPULATION_LOG    +
                  URBANIZACAO                    + EVANGELICOS       +
                  NAOBRANCOS_2010                + POBRES            + 
                  GINI                           + 
                  factor(REGIAO)                 +
                  BOLSONARO_1T_VOTOS_VALIDOS, 
                data=dados)

models$m3 <- lm(VOTOS_SEGURANCA_PROPORCAO ~ 
                  HOMICIDIO_2018_LOG*factor(CAPACIDADE)              +
                  CANDIDATOS_SEGURANCA_PROPORCAO + NEP               +
                  MAGNITUDE                      +
                  PREFEITOS_ESQUERDA             + PREFEITOS_CENTRO  +                  
                  ESQUERDA_2014                  + CENTRO_2014       + 
                  PIB_PERCAPITA_LOG              + POPULATION_LOG    +
                  URBANIZACAO                    + EVANGELICOS       +
                  NAOBRANCOS_2010                + POBRES            + 
                  GINI                           + 
                  factor(REGIAO), 
                data=dados)

models$m4 <- lm(VOTOS_SEGURANCA_PROPORCAO ~ 
                  HOMICIDIO_2018_LOG*factor(CAPACIDADE)              +
                  CANDIDATOS_SEGURANCA_PROPORCAO + NEP               +
                  MAGNITUDE                      +
                  PREFEITOS_ESQUERDA             + PREFEITOS_CENTRO  +                  
                  ESQUERDA_2014                  + CENTRO_2014       + 
                  PIB_PERCAPITA_LOG              + POPULATION_LOG    +
                  URBANIZACAO                    + EVANGELICOS       +
                  NAOBRANCOS_2010                + POBRES            + 
                  GINI                           + 
                  factor(REGIAO)                 +
                  BOLSONARO_1T_VOTOS_VALIDOS,
                data=dados)

# Command to create the table 1 of the article with the OLS models

# The results from the table bellow generated by this R command were incorporated 
# into Table 1 of the article, with only the variable labels adjusted 
# in the article (refer to the code dictionary for the associated labels).

stargazer(models, 
          out="./tabelas/table_1_ols_models.html", 
          type="html",
          title = "Modelos Aditivos e Interativos",
          dep.var.labels = "% de votos municipais a candidatos à Câmara Federal da área de segurança",
          column.labels = c("Modelo 1",
                            "Modelo 2",
                            "Modelo 3",
                            "Modelo 4"),
          align = TRUE )

####################################################
## 05. OLS Models - Cluster (UF)                  ## 
##                                                ##
## Article: Table A-2 (appendix)                  ##
####################################################

models_cluster <- list()

models_cluster$m1 <- feols(VOTOS_SEGURANCA_PROPORCAO ~
                             HOMICIDIO_2018_LOG             + factor(CAPACIDADE)  +
                             CANDIDATOS_SEGURANCA_PROPORCAO + NEP               +
                             MAGNITUDE                      +
                             PREFEITOS_ESQUERDA             + PREFEITOS_CENTRO  +                  
                             ESQUERDA_2014                  + CENTRO_2014       + 
                             PIB_PERCAPITA_LOG              + POPULATION_LOG    +
                             URBANIZACAO                    + EVANGELICOS       +
                             NAOBRANCOS_2010                + POBRES            + 
                             GINI                           + 
                             factor(REGIAO), 
                           data=dados)

models_cluster$m2 <- feols(VOTOS_SEGURANCA_PROPORCAO ~
                             HOMICIDIO_2018_LOG             + factor(CAPACIDADE)  +
                             CANDIDATOS_SEGURANCA_PROPORCAO + NEP               +
                             MAGNITUDE                      +
                             PREFEITOS_ESQUERDA             + PREFEITOS_CENTRO  +                  
                             ESQUERDA_2014                  + CENTRO_2014       + 
                             PIB_PERCAPITA_LOG              + POPULATION_LOG    +
                             URBANIZACAO                    + EVANGELICOS       +
                             NAOBRANCOS_2010                + POBRES            + 
                             GINI                           + 
                             factor(REGIAO)                 +
                             BOLSONARO_1T_VOTOS_VALIDOS, 
                           data=dados)

models_cluster$m3 <- feols(VOTOS_SEGURANCA_PROPORCAO ~ 
                             HOMICIDIO_2018_LOG*factor(CAPACIDADE)              +
                             CANDIDATOS_SEGURANCA_PROPORCAO + NEP               +
                             MAGNITUDE                      +
                             PREFEITOS_ESQUERDA             + PREFEITOS_CENTRO  +                  
                             ESQUERDA_2014                  + CENTRO_2014       + 
                             PIB_PERCAPITA_LOG              + POPULATION_LOG    +
                             URBANIZACAO                    + EVANGELICOS       +
                             NAOBRANCOS_2010                + POBRES            + 
                             GINI                           + 
                             factor(REGIAO),
                           data=dados)

models_cluster$m4 <- feols(VOTOS_SEGURANCA_PROPORCAO ~ 
                             HOMICIDIO_2018_LOG*factor(CAPACIDADE)        +
                             CANDIDATOS_SEGURANCA_PROPORCAO + NEP               +
                             MAGNITUDE                      +
                             PREFEITOS_ESQUERDA             + PREFEITOS_CENTRO  +                  
                             ESQUERDA_2014                  + CENTRO_2014       + 
                             PIB_PERCAPITA_LOG              + POPULATION_LOG    +
                             URBANIZACAO                    + EVANGELICOS       +
                             NAOBRANCOS_2010                + POBRES            + 
                             GINI                           + 
                             factor(REGIAO)                 +
                             BOLSONARO_1T_VOTOS_VALIDOS,
                           data=dados)

# Clusterizando em one-way (one-way clustering)
controle_cluster <-etable(models_cluster$m1,
                          models_cluster$m2,
                          models_cluster$m3,
                          models_cluster$m4, 
                          se = "cluster", cluster =c("UF"))

# Command to create the table A2 of the article with the OLS models

# The results from the table bellow generated by this R command were incorporated 
# into Table A-2 of the article, with only the variable labels adjusted 
# in the article (refer to the code dictionary for the associated labels).

write_xlsx(controle_cluster, "./tabelas/table_a2_ols_models_cluster.xlsx")

####################################################
## 06. Marginal Effect                            ## 
##                                                ##
## Article: Table 2                               ##
####################################################
library(margins)
library(broom)

marginal_effects <- margins(models$m3, 
                            variables = "HOMICIDIO_2018_LOG",
                            at = list(CAPACIDADE = c(0, 1))
)

# Convert to a tabular table with broom
marginal_effects_summary <- summary(marginal_effects)
print(marginal_effects_summary)

# Save the table as CSV file

# The results from the table bellow generated by this R command were incorporated 
# into Table 2 of the article, with only the variable labels and columns names adjusted 
# in the article (refer to the code dictionary for the associated labels).

write.csv(marginal_effects_summary, "./tabelas/table_2.csv", row.names = FALSE)

####################################################
## 07. VIF testls                                 ## 
####################################################

vifs <- list()

vifs$vif_m1 <- vif(models$m1)
vifs$vif_m2 <- vif(models$m2)
vifs$vif_m3 <- vif(models$m3,  type = 'predictor')
vifs$vif_m4 <- vif(models$m4,  type = 'predictor')

# Print VIFs
print(vifs)

# Save the VIFs - Table (models - m3 and m4)
stargazer(vifs, out="./tabelas/vifs_table.html", type="html")

####################################################
## 08. Descriptive Statistics                     ## 
##                                                ##
## Article: Table A-1                             ##
####################################################

## Article: Table A-1

# Create a subset of the dataframe
descritive_dados <- dados %>%
  dplyr::select("VOTOS_SEGURANCA_PROPORCAO","HOMICIDIO_2018_LOG","CAPACIDADE",
         "CANDIDATOS_SEGURANCA_PROPORCAO","NEP","MAGNITUDE","PREFEITOS_ESQUERDA",
         "PREFEITOS_CENTRO","ESQUERDA_2014","CENTRO_2014","PIB_PERCAPITA_LOG",
         "POPULATION_LOG","URBANIZACAO","EVANGELICOS","NAOBRANCOS_2010","POBRES","GINI")

# Use the command bellow to measure the descriptive statistics

# The results from the table below generated by the `skim` command were incorporated
# into Table A-1 of the article, which was generated directly in the article's text file.
# The variable labels in the article were adjusted according to 
# the codebook for the associated labels used in this table.

# (Table A-1) Descriptive statistics
table_a1 <- skim(descritive_dados)
print(table_a1)
write.csv(table_a1, "./tabelas/table_a1.csv", row.names = FALSE)

####################################################
## 09. Spatial Data Visualization                 ## 
##                                                ##
## Article: Figure 1                              ##
## Article: Figure A-1 (appendix)                 ##
####################################################

# Import the municipalities shapefile to R
shp_cities <- read_sf("./shapefiles/BR_Municipios_2022.shp")
class(shp_cities) # To check the file
shp_cities        # To check the file

# Merge the shapefile with the dataset
dados$CD_MUN    <- as.character(dados$COD_MUN_IBGE_SETE)
dados_map       <- left_join(shp_cities, dados, by = "CD_MUN")

# Remove "Fernando De Noronha"
dados_map <- dados_map %>% filter(CD_MUN != 2605459)

## Article: Figure 1
# (Figure 1) Plot the dependent variable
figure_1 <- ggplot(dados_map) +
  geom_sf(aes(fill = VOTOS_SEGURANCA_PROPORCAO),
          color = "transparent") +
  scale_fill_gradientn(colors = c("white", "lightcoral", "red", "darkred"),
                       values = scales::rescale(c(0, 0.25, 0.5, 1))) + 
  labs(fill = "%") +
  theme_void() +  
  theme(legend.key.width = unit(0.5, "cm"),
        legend.position = "right")  
# Te save the figure 1
ggsave("./figuras/figure_1.png", plot = figure_1, width = 10, height = 6, dpi = 300)

# Article: Figure A-1 (appendix)
# (Figure A-1) Plot the variable "Proporçao de Candidatos da Seguranca"
figure_a1 <- ggplot(dados_map) +
  geom_sf(aes(fill = CANDIDATOS_SEGURANCA_PROPORCAO),
          color = "transparent") +
  scale_fill_gradientn(colors = c("white", "lightcoral", "red", "darkred"),
                       values = scales::rescale(c(0, 0.25, 0.5, 1))) + 
  labs(fill = "%") +
  theme_void() +  
  theme(legend.key.width = unit(0.5, "cm"),
        legend.position = "right")  
# Te save the figure 2
ggsave("./figuras/figure_a1.png", plot = figure_a1, width = 10, height = 6, dpi = 300)
